Unveiling environmental entanglement in strongly dissipative qubits 
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The coupling of a qubit to a macroscopic reservoir plays a fundamental role in understanding 
the complex transition from the quantum to the classical world. Considering a harmonic environ- 
ment, we use both intuitive arguments and numerical many-body quantum tomography to study 
the structure of the complete wavefunction arising in the strong-coupling regime, reached for in- 
tense qubit-environment interaction. The resulting strongly-correlated many-body ground state is 
built from quantum superpositions of adiabatic (polaron-like) and non-adiabatic (antipolaron-like) 
contributions from the bath of quantum oscillators. The emerging Schrodinger cat environmental 
wavefunctions can be described quantitatively via simple variational coherent states. In contrast to 
qubit-environment entanglement, we show that non-classicality and entanglement among the modes 
in the reservoir are crucial for the stabilization of qubit superpositions in regimes where standard 
theories predict an effectively classical spin. 



The study of dissipative quantum phenomena, namely 
the interaction of a quantum object (a qubit) with an in- 
finite number of environmental degrees of freedom, lies at 
the frontier of modern science and technology, with deep 
implications for fundamental quantum physics^, quan- 
tum computing^, and even biolog}'^^. While quantum 
information stored in the qubit subsystem is lost during 
the coupling with the unobserved degrees of freedom in 
the reservoir, it is in principle preserved in the entan- 
gled many-body state of the global system. The precise 
nature of this complete wavefunction has received little 
attention, especially regarding the entanglement gener- 
ated among the reservoir states. Our purpose here is to 
unveil a simple emerging structure of the wavefunctions 
in open quantum systems, using a complementary com- 
bination of numerical many-body quantum tomography 
and a novel analytical variational theory. 

An archetype for quantitatively exploring the quantum 
dissipation problenP^ is to start with the simplest quan- 
tum object, a two- level system describing a generic quan- 
tum bit embodied by spin states {| t)^ I i)}? and to couple 
it to an environment consisting of an infinite collection of 
quantum oscillators a\ (with continuous quantum num- 
ber k and energy hwk)- Quantum superposition of the 
two qubit states is achieved through a splitting A act- 
ing on the transverse spin component, while dissipation 
(energy exchange with the bosonic environment) and de- 
coherence are provided by a longitudinal interaction term 
Qk with each displacement field in the bath. This leads to 
the Hamiltonian of the celebrated continuum spin-boson 
model (SBM)^^: 

H = - X! y (4 + + (1) 
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where we set = 1, and the sums can be considered as 
integrals by introducing the spectral function of the envi- 
ronment, J(co') = X]/c^fc^(^ ~ ^k)' The generality of the 
SBM makes it a key model for studying non-equilibrium 
dynamics, non-Markovian quantum evolution, biological 
energy transport, and the preparation and control of ex- 
otic quantum states in a diverse array of physical and 
chemical systems^ ^. 

The possibility of maintaining robust spin superposi- 
tions in the ground and steady states of the SBM has 
attracted considerable attention, prim arily due to its im- 
plications for quantum computin^^^. Previous numeri- 
cal approaches have hitherto mainly focused on observ- 
ables related to the qubit degrees of freedonPHm^ whilst 
a description of the global system-environment wavefunc- 
tion has been confined to simpler variational studies"*^^-'^. 
This variational theory readily predicts the formation of 
semiclassical polaron states, which involve the adiabatic 
response of the environmental modes to the spin tun- 
neling. Strong entanglement between the qubit and the 
bath is generated in this process. We shall demonstrate 
here that the many-body ground state of Hamiltonian ([T]) 
contains additional non-classical correlations among the 
environmental oscillator modes arising from their non- 
adiabatic response to the spin-flip processes. These new 
non-classical contributions to the wavefunction are key 
for the actual stabilization of qubit superpositions rel- 
ative to the semiclassical picture, and naturally emerge 
from a variational framework beyond the adiabatic po- 
laron approximation. 

In order to enlighten the nature of these emergent 
non-classical environmental states, we first analyse the 
SBM by performing the (unitary) polaron transforma- 
tion H = U HW , where U exp{ -a^ ^ (a^ - ) } , 
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which removes the hnear interaction term in Eq. Q. 
This transforms the Hamiltonian to a basis in which oscil- 
lator wavefunctions are displaced according to the 2:-axis 
projection of the spin: 



H 



-h.c. - 
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where Eji = g'^/ (4co'/e) is the reorganisation energy of 
the bath. For A = 0, the ground state of H is doubly 
degenerate, and is given by the product of the bosonic 
vacuum and the spin states, |^^,o) = | t) ^ |0) and 
|^i,o) = |i) |0), in the transformed basis (denoted by 
tildes). It thus corresponds to polaronic wavefunctions 
in the original frame, where the positive/negative sign of 
the displacement is fully correlated to the spin projec- 
tion {adiabatic response): |^t,^fc/2u;fc) = 1 1) I -^Qk/'^^k) 
and |^;,-^,/2u;J = I ;) (8) I - Qk/'^^k)- The two-fold de- 
generate ground state thus takes the form of a prod- 
uct of semiclassical coherent states (displaced oscilla- 
tors) I ± /fc) = e'^^kfk{a^-a^)^^\^^ ^-^^Yi displacements 
fk = ±^/c/2a;/e which shift each oscillator to the min- 
imum of its static spin-dependent potential. This po- 
tential is evident in Eq. ([T]) for A = and is shown 
explicitly in Fig. [T]A.. In the presence of spin tunneling 
(A 7^ 0), one needs to understand the effect of the op- 
erators K± = Aa^e^^ki9k/uJk)ial-a,^) ^ ^j^-^j^ 

correlate spin flip processes with simultaneous displace- 
ments of all oscillator states. As we now show, these 
correlations ultimately control the ground state qubit su- 
perposition. 

Polarons, antipolarons, and ground state ansatz. The 
optimum oscillator displacements result from a competi- 
tion between the two terms appearing in Hamiltonian ([2| , 
namely spin tunneling A versus oscillator kinetic energy. 
Within the transformed frame, consider the coupling in- 
duced by the tunneling operator K± between one of the 
doubly degenerate states, say |^|,o) = I i) (8) |0), and 
a spin-flipped state with arbitrary displacement function 
Ik = fk-gk/^^k^ l^tJfc^ ^ I t)^\fk) {fk is the displace- 
ment in the original frame). The matrix element for this 
is 



(3) 



The elastic displacement energy of oscillator /c in | ) 

is {'^^jjukalakl'^^jj = ujkfl- The polaron trans- 
formed Hamiltonian reveals the inherent competition be- 
tween the (elastic) energetic cost of mixing displaced os- 
cillators into the ground state, favouring = 0, and 
the exponential suppression of the spin kinetic energy, 
given by the reduced tunneling matrix element in Eq. ([s]), 
which rather favours fk = —Qk/^k- For high-frequency 
modes, the elastic energy cost dominates and tunneling 
between spin states is governed by environment states 
with //e = 0, gaining only a small (renormalized) tun- 
neling energy A^^ = /S.e~^^ki9k/^k) ^ fQ^. strong 



qub it-environment interaction. In the original frame, the 
corresponding displacement is fk = -^Qk/'^oJk^ which im- 
plies that these 'fast' oscillators instantaneously (adia- 
batically) tunnel with the spin between the minima of 
their elastic potentials - see Fig. [l] In the opposite limit 
of low- frequency modes {uok <C A^^), the elastic energy 
barrier is weak; mixing between spin states is instead 
governed by the matrix element ([3|. Returning to the 
original frame, one gets an energy gain of the hare tunnel- 
ing energy A when fk = —ghj^uj^. As shown in Fig.[l]C, 
this corresponds to spin tunneling with non- adiabatic re- 
sponse of the oscillators, which are displaced in the op- 
posite direction from the adiabatic modes. We naturally 
dub these contributions to the wavefunction antipolaron 
states. At intermediate frequencies, we expect that both 
polaronic and antipolaronic responses occur, leading to a 
two-polaron ansatz for the ground state (in the original 
frame): 



(4) 
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with p the relative weight of the polaron and antipo- 
laron components. Note that this ansatz fully respects 
the symmetries of the Hamiltonian. 

This state reduces to standard (adiabatic) polaron the- 
ory when p = and /^°^' = gk/2ujk^ and to the varia- 
tional polaron state of Silbey and Harris (SH)?^!^ when 
p = and the function /^°^" is varied to minimise the to- 
tal ground state energy E = (GS'^p^^- IG^^p^^-). As we 
shall compare our ansatz Q to these simpler theories, a 
brief description of them is given in the Supplementary 
Information. For p 7^ 0, the environment wave function 
for each spin projection is a multi-modal Schrodinger cat 
state involving a superposition of polaronic and antipo- 
laronic components, leading to considerable mode entan- 
glement. The critical observation is that such superposi- 
tion of displaced states lowers the energy of the ground 
state by stabilising the spin energy. For the state Q the 
spin tunneling energy Et = (A/2)(G5'2p°^- |cr^|G!s^P°^-) 
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The first two terms reflect an exponentially suppressed 
renormalized tunneling rate. Indeed, for strong coupling, 
the displacements /^°^' and /^"*^" are large, and the as- 
sociated contribution to the spin energy becomes vanish- 
ingly small. However, the overlap between the polaron 
and antipolaron contributions (the third term) will not 
be suppressed if, as we expect, /^^ 



-/r'. Thede- 



velopment of a small but finite antipolaron weight p thus 
allows the environment to minimise its displacement en- 
ergy whilst maintaining significant overlap between the 
environment-dressed spin states. 

Single mode. Before tackling the challenging many- 
mode situation, we develop intuition about the polaron- 
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FIG. 1. Origins of polaron and antipolaron displacements in environment wavefunctions. In plots A-C, black 
dashed lines are the spin-dependent potential energies of a single harmonic oscillator in the absence of spin tunneling [see 
Eq. ([T])], while blue (red) curves are the gaussian wavefunctions (in real space x) of the oscillator on the (cr^) = 1 (—1) potential 
surfaces. A. Polarons. For a high frequency mode {uj ^ A), transitions to other oscillator states on the potential surfaces 
are suppressed by the steep curvature of the potentials; oscillator displacement adiabatically tunnels with the spin between 
minima of the potentials, suppressing the tunneling amplitude by the reduced overlap of the displaced oscillator wave functions 
to a value A^. B. Non-adiabatic response in Silbey-Harris variational polaron theory. Low frequency modes 
(cj ^ A) have shallow potentials, leading to well-separated minima. Poor wave function overlap prevents tunneling of the spin 
between minima, destroying spin superposition. Variationally-determined displacements adjust to smaller values, sacrificing 
their displacement energy to maintain the spin-tunneling energy through better overlap. C. Antipolaron response of non- 
adiabatic oscillators. For modes with cj ~ A, spin flips that do not change the position of the oscillators (and thus have 
unsuppressed amplitude A) may become low enough in energy to compete with the overlap-suppressed inter-minima tunneling. 
The oscillator wavefunctions correlated with spin are now superpositions of displaced coherent states with opposite signs. D. 
Ground state wavefunction components of a single oscillator. Spin up (blue) and spin down (red) components are 
shown for the exact ground state (circles), our variational polaron-antipolaron state (solid lines), and the Silbey-Harris ansatz 
(dashed lines) [A/ui = 4, g/oJi = 3]. The exact result shows distinct antipolaron features which are well captured by the 
variational polaron-antipolaron state. The Silbey-Harris ansatz shows reduced displacements and thus poor agreement with 
the exact result. 



antipolaron ansatz in the simplest case of a single envi- 
ronmental mode with energy ui and coupling gi. This 
case is easily diagonalised numerically (see also Ref. [24] 
for an exact solution); note that a similar ansatz for the 
single- mode Rabi model (without reference to pola ron 
theory) has been previously explored numeric alljP^^. In 
Figure [ijl) we compare the spatial wavefunctions of the 
oscillator correlated with each spin state with those ob- 
tained from the ansatz Eq. Q following a numerical op- 
timization of p, /f°^', and f^^^^- to minimise the ground 
state energy. Choosing oscillator parameters where we 
expect non-adiabatic response, namely uji < A, we find 
that both wavefunctions clearly show a superposition of 
polaron and antipolaron contributions, with much larger 
displacements compared to the prediction of the SH the- 
ory (single polaron case p = 0). The agreement of the 
diagonalised and the two-polaron ansatz ground state 
wavefunctions is extremely good, as well as the energies 
and spin observables, even for a coupling strength as large 
as ^ = 3ui (see also Supplementary Information). As 
motivated above, the emergence of an antipolaron com- 
ponent in the environment enhances the overlap of the 
tunneling states. The single polaron SH state fails in 
this regard (see Figures [l^ and ID, and Supplementary 
Information) as it finds itself frustrated between mini- 
mizing the elastic energy and maintaining good overlap 
between the opposite spin states: the resulting displace- 



ments are thus totally wrong. 

Two-mode antipolaron entanglement. Having con- 
firmed the emergence of non-adiabatic antipolaron con- 
tributions in the case of a single mode, we now consider 
the case of a two-mode SBM and, in particular, test our 
proposal Eq. Q that the two-mode wave function dress- 
ing a given spin state will be entangled. 

Fig. |2] shows the spin- up component of the two- mode 
wavefunction as a function of the two independent spa- 
tial coordinates of the modes (xi and X2) for two modes 
taken at different frequencies UJ2 = 2cJi. The ground 
state wavefunctions were determined by exact numerical 
diagonalisation. We see the clear development of an an- 
tipolaron component to the wavefunction (Fig. |2^) for 
low-energy non-adiabatic modes, in contrast to the sit- 
uation of high-energy adiabatic modes (Fig. [2]A.). How- 
ever, we see that only two peaks appear in the wave- 
function - those along the diagonal line Xi = X2 - indi- 
cating unambiguously that this two-mode wavefunction 
takes the inter- mode entangled form |/f°^") l/l^^') + 
^lyanti.^ (g) l/f"*' ). This cau be contrasted with a hy- 
pothetical polaron-antipolaron product state {|/f°^') + 

^lyanti.^l ^ {\fP^^-) + j ^^^^^ ^^^^^ ^^^^^^ ^-g, 

play four peaks, as shown in Fig. [2]C. The implications 
of this inter-mode entanglement for the entropy of the 
reservoir modes is given in Supplementary Information. 
Again, one can check that the variational energy of the 
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FIG. 2. Two mode wavefunctions. A-B. Contour plots in real space of the spin-up projected joint oscillator wavefunctions 
of two modes obtained from exact diagonalisation. A. Polarons. For high frequency modes (002 = 2uji = 0.04 > A = 0.01), the 
wavefunction is a single, displaced gaussian, in qualitative agreement with Silbey-Harris theory. B. Entangled antipolarons. 
Low frequency modes {002 = 2uji = 0.004 < A = 0.01) show the development of an antipolaron component, visible in the 
(Xi < 0^X2 < 0) quadrant, in addition to the Silbey-Harris state. C. Product state. Hypothetical wavefunction obtained 
from a product state of polaron-antipolaron superpositions for each mode, showing symmetric off-diagonal peaks. These features 
are absent in B, indicating that the exact joint wavefunction is not a product state but, in contrast, is entangled as described 
in the text. 



two-mode ground state is remarkably close to the exact 
energy. 

Multi-mode spin-boson model. We now turn to the 
more challenging many-mode situation, tackling the con- 
tinuum spin-boson model ([T]). A direct diagonalisation 
of the Hamiltonian is now hopeless; however, recent 
computational progress has opened the way to calcu- 
lating ground state averages of arbitrary operators, for 
instance using the bosonic Numerical Renormalization 
Group (NRG)^^, which we will use to test the general- 
ized polaron state Q . A key feature in the NRG method 
is the use of a logarithmic discretization of the energy 
spectrum of the bath, which ensures the stability and 
convergence of an iterative diagonalization of the impu- 
rity model^^. In order to directly compare with the vari- 
ational results, we use the same discretization in defining 
the polaronic ansatz Eq. Q, incorporating the changing 
measure in ujk into the definitions of the fk. 

We focus here on the standard case of ohmic 
dissipatiorl^, although our following results should ap- 
ply similarly to other types of spectral density. The con- 
tinuous bath of bosonic excitations assumes then a lin- 
ear spectrum in frequency, J{uj) = 2auj6{uJc — cj), up to 
a high energy cutoff ujc and with dimensionless dissipa- 
tion strength a. Weakly damped Rabi oscillations of the 
qubit for a <C 1 are known to completely fade away in 
the strong dissipation regime a > 0.4, where the qubit 
becomes strongly entangled with its environment. The 
bare qubit frequency A is heavily renormalized in this 
regime to the smaller value A^^ = A{Ae/uJc)'^^''^~^\ for 
A/uJc <C 1, which can thus be driven to zero for the crit- 
ical dissipation strength o^c — 1, indicating a quantum 
critical point. 



As a first step towards understanding the many-mode 
situation, we consider the variational solution obtained 
from the two-polaron ansatz Q (the variational equa- 
tions are given in the Supplementary Information). This 
leads to the polaronic and antipolaronic displacements 
shown in Fig. [3]A, which exemplify the physical picture 
introduced above (see especially Fig.[T]): polaron and an- 
tipolaron states show equal and opposite displacements 
at low energies (typically for Uk <C cjc), but merge to- 
gether to produce a fully polaronic state at high en- 
ergy, where the environment responds adiabatically to 
the spin. The variational theory is thus able, without ad- 
ditional physical input, to generate the correct crossover 
from non-adiabatic to adiabatic behavior of the antipo- 
laron component with increasing energy. 

The presence of the antipolaron component has a large 
impact on the ground state spin average: Fig. [3^ com- 
pares the result of the one- and two-polaron variational 
states to that computed with NRG (numerically exact 
result, used as a benchmark). In the one-polaron (SH) 
limit {p = 0), one finds readily —{(Tx) = Ar/A = 
{Ae/uJc)'^^^^~^\ which incorre ctly vanishes at the crit- 
ical dissipation strength ac = j "*^^* I On the other hand, 
the emergence of antipolaron correlations at low energy, 
namely /^"*^" — — /^°^" for ^ uJc, helps in maintaining 
a finite value for {ax), due to the perfect cancellation of 
the displacements within the exponential in the last term 
of Eq. ([5|. This success of the antipolaron ansatz Q is 
illustrated in Fig. [3^. 

Our objective now is to demonstrate the peculiar 
inter-mode entanglement properties of the antipolaron 
ansatz 0. While one cannot plot the complete many- 
body wavefunction in the case of many environmental 
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FIG. 3. Displacements and spin average in the many-mode case. A. Displacements determined variationally from 
the two-polaron ansatz Eq. Q, showing the emergence of an antipolaron component for low energies, with equal and opposite 
displacement to the polaron state. The antipolaron state merges smoothly onto the polaron state at high energy as the 
adiabaticity of the oscillators with respect to tunneling of the spin is recovered (the NRG logarithmic discretisation of the bath 
spectrum is used here, namely frequency points are evenly spaced on a logarithmic scale; Note that a point at higher energy 
is associated to a larger energy window of the continuum spectrum, leading to the saturation of /^°^' for high frequencies, 
instead of the fall off obtained for a linear energy mesh). [Parameters: a — 0.5 and A = 0.01.] B. Ground state averaged spin 
amplitude —{ctx) as a function of dissipation strength a computed with the NRG (circles) for A/ujc = 0.01, and compared to 
the one-polaron (red line) and two-polaron (blue line) predictions. A clear breakdown of the one-polaron Silbey-Harris ansatz 
occurs at strong dissipation, while the two-polaron trial state accounts for the correct behavior up to the quantum critical point 
(ac = 1), due to preserved tunneling amplitude via the antipolaron component of the wavefunction. 



modes, a useful strategy to assess the validity of the 
trial state Q lies in recent interest in quantum tomogra- 
pj^^^mmS^ wherein the reduced density matrix in a smaller 
projected Hilbert space is fully characterized. For the 
problem at hand, we trace out all modes except the qubit 
degree of freedom together with an arbitrary bath mode 
with given quantum number k] this defines a spin and k- 
mode excluded environment denoted "env/spin-h/c" . The 
reduced ground state density matrix in the joint qubit 
and /c-mode subspace reads 



in the regime of strong dissipation (a > 0.5): 



Pspin+/c — Trenv/spin+/c|G5')(GS'|. 



(6) 



We focus here on the off-diagonal part (with respect to 
the qubit axis of quantization) of the Wigner distribu- 
tion associated to this density matrix as a function of the 
classical displacement X. We expect on physical grounds 
that this component will be most sensitive to the antipo- 
laronic correlations. Its standard definition i^ 



^Cr^Pspin+fc 



(7) 

see Methods for the NRG implementation and Supple- 
mentary Information for discussion of the spin-diagonal 
part of the Wigner distribution, which emphasizes in- 
stead the polaronic part of the total wavefunction. From 
the two-polaron trial state Q and equation ([7|), it is 
straightforward to find the form of the Wigner function 



1 / ppol. I ranti. \ 



(8) 



^pol. panti. 



pol. 



For high energy modes uj^ ^ ^c-, 

W^^^{X) should show 
a single peak centered around X = 0, as both polarons 
adiabatically follow the spin tunneling so that /^"*^" be- 
comes close to the polaron displacement /^°^' . For modes 
of lower energy, antipolaron displacements emerge, and 
the peak separates into two lobes with displacements 
±[/r^" - - These simple predictions of 

the two-polaron variational state are clearly seen in the 
numerical NRG result in Fig. [4j strongly supporting the 
existence of non-adiabatic oscillator states in the envi- 
ronment. 

We finally wish to assess more directly the entangle- 
ment among the environmental states that is suggested 
by the antipolaron ansatz Eq. (4). For this purpose, we 
consider the entanglement entropy S'spin+k of the joint 
spin and /c-mode subsystem with respect to the other 
modes of the bath: 

^spin+k = — Trspin+/c [Pspin+/c log Pspin+/c] 7 (9) 

with the reduced density matrix defined in Eq. (|6|. We 
also introduce the spin entanglement entropy 6'spin = 
-Trspin [Pspinlogpspin] with Pspin = Tr^ [pspin+/c] • The 
difference between these two quantities can be computed 
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FIG. 4. Quantum tomography in the many-mode case. 

Transverse Wigner distribution defined in Eq. ([7|, as ob- 
tained from NRG, for two different modes, one at fiigfi energy 
cjfc ^ Ai^ (top red curve) and tfie otfier at intermediate en- 
ergy cjfc ^ Ai^ (bottom blue curve). A decomposition of the 
intermediate energy case into two shifted Gaussians is per- 
formed (thin blue lines) according to Eq. ([s]), showing the 
emergence of antipolaron correlations. [Parameters a = 0.8 
and A/cjc = 0.01.] 



from the NRG (see Methods) and is plotted as a func- 
tion of mode frequency in Fig. |5| For small dissipation, 
a < 0.5, this entropy difference is mostly negative, as ex- 
pected from the correlations built into the Silbey-Harris 
state, which consist only of non-entangled environmen- 
tal states within each spin-projected component of the 
wavefunction (see Supplementary Information). In con- 
trast, at strong dissipation, this entropy difference be- 
comes positive and shows a strikingly large enhancement 
near the scale A^^. The excess entanglement entropy, 
above that of the spin alone, comes from entanglement 
within the bath of oscillators. This is a sensitive signa- 
ture, then, that the spin projected wavefunction is not 
simply a product of oscillator states as in the SH ansatz 
but rather involves substantial entanglement. The na- 
ture of this entanglement in the simpler two-mode case 
is explored further in the Supplementary Information. 
Note especially the large energy window where the en- 
tropy peak develops: the excess entanglement spreads 
from low to high frequency modes due to the massive 
entangling power of the spin tunneling operator K-^ dis- 
cussed above in the polaron basis H of Eq. ([I]). The 
existence of inter-mode bosonic correlations on a wide 
energy range makes also an interesting connection to the 
underlying (although hidden in the spin-boson model) 
fermionic Kondo physic^^^. 

In conclusion, we have shown how antipolaron contri- 
butions emerge in the ground state wavefunction of the 
spin-boson model, causing non-classical Schrodinger cat- 
like environmental states. The approach here can be used 
as a general framework to expand and rationalize many- 
body wavefunctions in strongly interacting open quan- 
tum systems. Experimentally, proposals to realize the 
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FIG. 5. Joint entanglement entropy in the many-mode 
case. The joint entanglement entropy of the qubit and a 
given (jjk mode, defined by Eq. Q, is calculated with NRG for 
A/cJc = 0.01. The entropy of the qubit alone is subtracted. 
Negative correlations for a < 0.5 are related to the Silbey- 
Harris one-polaron state, while the strong positive peak can 
only be accounted for by entanglement among the modes of 
the bath. This excess entropy is, then, a sensitive measure 
of the subtle non-classical correlations among the bath modes 
that are generated by the coupling to the qubit. 



strongly dissipative spin boson model using a supercon- 
ducting qubit coupled to arra^ of Josephson junctions 
have been made very recently^^^*^. The recent pro gress 
in quantum tomography of superconducting qubit 9^^^ 
raises thus the challenge to measure in such setups the 
massive entanglement of environment oscillators that was 
unveiled here. The present work offers several direc- 
tions for future research, especially the generalization to 
time-dependent phenomena such as the study of quan- 
tum quenches and spin dynamics at strong dissipation, 
where standard (weak-coupling) Bloch-Redfield theory^ 
is known to fail. 

Methods 

The numerical solution of the few-mode spin-boson 
Hamiltonian ([T]) relies on standard diagonalisation pro- 
cedures. In the case of a continuous bath of oscilla- 
tors, a different strategy is used. First, logarithmic 
shell blocking of the bosonic modes onto energy inter- 
vals [A~'^~'^ujc, A~^uJc\ with A = 2 is performed: 



4 = / dkal. 



(10) 



The resulting discrete Hamiltonian, which spans from ar- 
bitrarily small energy up to the high energy cutoff ujc, is 
then iteratively diagonalised according to t he Nu meri- 
cal Renormalization Group (NRG) algorithnjm^. The 
novel part of the simulations performed for this work lies 
in the computation of the Wigner distribution reduced 
to the joint spin and single /c-mode subspace. In order to 
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implement Eqs. (|6]|7|, we first define arbitrary moments 
of the chosen oscillator k of frequency tOk = uJc-^~^: 



A 



(GS\aAair\a„ 



(11) 



with i = O^x^y^z labelling the Pauli matrices related to 
the spin projection (we take ctq = 1) and m, posi- 
tive integers. Such ground state observables are read- 
ily computed within the NRG algorithm (for typically 
< m, < 10). One can then expand Eq. ([7|) in a 
power series in A and A, yielding 



/ V (7i\m^m' 



m\m'\ 



(12) 



The Wigner distribution is now solely expressed in terms 



of the NRG-computable moments 

A similar strategy is used for the computation of the 
entanglement entropy ([9| from the reduced density ma- 
trix />spin+/c, which acts within the subspace spanned 
by the qubit and a single bosonic mode k. We start 
by defining the joint spin and Fock projection operator 

^^a-mm' ~ l^/c ) (^fc | , SO that matrix elements of the 
ground state density matrix simply read 



(13) 



This quantity is a ground state average, hence readily 
computable by letting the operator O^..^ ^, evolve along 
the complete NRG flow. The eigenvalues of the matrix 
P^a^m m' ^^^ow ouc, finally, to obtain the desired entan- 
glement entropy. 

A last new piece of our work is the multipolaron gen- 
eralization of the previous single-polaron trial statd^^^ 
this is key for capturing easily the emergent non- 
adiabatic physics at strong dissipation. The variational 
method is straightforwardly implemented in the few- 
mode case by minimizing the average Hamiltonian ([T]) 
while using the double-polaron ansatz (4|. In the many 
mode case, despite having two sets of unknown functions, 
/^°^" and /^"*^', labeled by the continuous momentum /c, 
one can show that their form as a function of k is uniquely 
fixed from the variational principle, leaving a finite set of 
effective parameters to be determined (see Supplemen- 
tary Information). One finds that the displacement as- 
sociated with the first polaron follows qualitatively the 
standard behavior /^°^" = 0.5gk/{uJk + A^^) known from 
SH theor}^^^, with some quantitative deviations due to 
the feedback of the antipolaronic state /^"*^". The latter 
takes the approximate form /^"*^" — /^^^ -^^^ with a 
new energy scale Q that controls the crossover from non- 
adiabatic to adiabatic behavior as a function of mode 
energy (see Supplementary Information for the complete 
expression). The antipolaronic (non-adiabatic) charac- 
ter at low energy of the second contribution in the trial 
state Q is thus automatically guaranteed by the varia- 
tional principle. 
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Supplementary information for 
"Unveiling environmental entanglement in strongly dissipative qubits" 



We present here additional technical details and extra results, supporting the multi-polaronic description of the 
many-body ground state of the spin-boson model at strong coupling. We first consider the general two-polaron 
variational formalism for an arbitrary number of modes. The ground state energy and wavefunctions are then 
investigated for a wide range of parameter in the single-mode (Rabi) model, highlighting the emergence of an- 
tipolaron correlations and the possible breakdown of the single-polaron Silbey-Harris ansatz. The two-mode Rabi 
model is afterwards considered, with the emphasis on entropic issues, which provide interesting signatures of envi- 
ronmental entanglement. The need for additional antipolaronic contributions in the wavef unction is also discussed. 
Finally, the continuous spin-boson model is further explored, with detailed derivations of the Wigner functions 
pertaining to the reduced qubit and single mode Hilbert space, as well as extra comparisons between Numerical 
Renormalization Group simulations and the variational technique. 

I. GENERAL TWO-POLARON VARIATIONAL FORMALISM 

A. Energetics 

We consider the unbiased spin-boson model^ ^, as defined by the Hamiltonian (1) of the main text: 

H = ^cr^ + ^a;/e4a/e-cr^^^(4 + a/e), (1) 
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with tunneling energy A, a set of oscillator frequencies cj^, and system-oscillator coupling strengths (assumed 
real). Here, cTz = |t) (t| ~ \\) Ul? with spin basis states \\) and and a\ (a/c) is the oscillator creation (annihilation) 
operator for mode k. Hamiltonian ([T]) spans the cases of few discrete modes up to a continuum of bosonic fields, in 
which case the discrete /c-sum ought to be replaced by an integral over energy. 
Our two-polaron variational ground state ansatz takes the form 



\Gs^^°'-) = It) [pi +fr') +P2 i+/r')j - \i) [pi -/r'-) +P2 i-/r')j , (2) 

where the bosonic part of the wavefunction involves coherent states of the form 

\±f,)=e^^^f^^-i--'^\0), (3) 

defined as products of displaced states, where |0) represents all oscillators being in the vacuum state. The presence of 
a Z2 symmetry, namely (It)^!!)^ \ i) ^ \ t )^ (^k ^ — <^/e), and the need for minimizing the spin tunneling energy 
enforces the chosen relative sign between the up and down components of the ground state wavefunction in Eq. ([2|. 
Both functions f^""^' and /^"*^- are taken as free parameters, and will be varied to minimise the total ground state 
energy E = (GS'^p^^' | H \ GS'^^'^^'). In contrast to the usual Silbey-Harris state (for which p2 = 0)^ ^, this more fiexible 
ansatz allows for the possibility of a superposition of variationally determined displaced oscillator states associated 
with each spin projection. 

Normalisation of [G^^p^^ ) implies the condition 

2pl + 2pl + 4j9ip2e"^ ^^^-^^"''"-^^"'"^^ = 1, (4) 
while the variational ground state energy is given by 

k 

-2j2g, +pi/r- +PiP2(/r- +/r")e-^^^'^^^°' -^^^""■)^) . (5) 

k 

In the limit that P2 ^ (and so pi 1/a/2) we recover the Silbey-Harris variational ground state energy, 

Esn = -^e-'^^^f^^' +YMfrf - E^^^/r'-. (6) 

k k 

while further setting ak = gk/(2ujk) in Eq. (|6| gives the (non- variationally optimal) bare polaron ground state energy 

£;poL = -^e-i^^«^/-^-E/^- (7) 

k 

Going back to the two-polaron variational state of Eq. we find that the ground state coherence is given by 

(o-^) = -2 (^p\e-'^^kUr'-? +p2g-2Efc(/r'")' -^2pip2e~^^^^^^^''''^^^'^'^^^ , (8) 

while the magnetisation (a^) = by symmetry in absence of magnetic field along (unless one enters the polarized 
phase at a > 1 in the ohmic spin-boson model). 



B. Variational displacements 

The two sets of displacements /^°^" and /^"*^" are variationally determined from the total energy E of Eq. ([s] 
according to dE/df^^^' = and dE/df^^^^- = 0, which gives the closed form: 

.poi. ^ 9k A.jpioj, + - A2[piP2(/r'i/r'-) + A12] 
" 2 (pjojk + Ai)(piwfc + A2) - biP2(/r'i/r') + Ai2]2 

.anti. ^ 9k AM^k + Ai) - ^l[piP2(/r'l/r"-) + A12] 

2 (p2^, + Ai)(piwfe + A2) - biP2(/r'i/r') + Ai2]2 ' 
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which is vahd for an arbitrary number of oscihator modes. Hence the generic /c-dependence of the displacement 
is fully constrained by the variational principle, which leaves a finite set of effective parameters to be determined 
self-consistent ly according to: 



Ai 


= Api{ - /r'l/r) + - /r'l/r') +PiP2(-^ +^)</ri/r') 


(11) 




— At)'^/ f 1 f \ -\- ^ T)-, TOO / f 1 f \ -L r)-,r)r.{ i~i hM f P^^' 1 f \ 

— ^P2\ Ik \Jk 2^^^^\ '"^^ ^^9)\Jk \Jk / 




A 


^ / i*pol. 1 ^anti. \ 1 _ _ /- ~\ / ^pol. 1 ^anti. \ 

= ^PiP2{-fk \Jk } +PiP2{(^ - g){fk \Jk } 


(13) 


Co 


=pi+'^PiP2{fr'vr'') 
=pi+2pMfr'vr') 

- Z^^kJk Ik 

k 


(14) 
(15) 
(16) 
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(17) 



k 



The one-polaron Silbey-Harris displacement /f^ = O.bgk/[(j0k + A( — /^'l/r')] trivially recovered from Eq. ([9| 
by letting P2 = 0. One can also check that /^"*^" — /^^ for /c ^ cx), while /^"*^" — — /^^ for /c ^ in the limit of 
strong dissipation. Thus the antipolaron displacement satisfies the expected adiabatic/non-adiabatic crossover as a 
function of energy cj/^, and this physical picture is naturally incorporated in the variational theory. 



II. SINGLE-MODE RABI MODEL: CHECKING THE WAVEFUNCTION 



It is illustrative to consider the simplified case of a single- mode within the environment, namely the Rabi model 
(see Ref. ^and references therein). In this situation, the model Hamiltonian may be diagonalised straightforwardly 
(numerically) and so the regimes of validity of our polaron- antipolaron ansatz, as well as of the Silbey-Harris and 
non-optimal bare polaron states, may be assessed. The Hamiltonian now becomes 



Hi 



-dx + 0Jia[ai 



3 ( t 



-ai). 



(18) 



with ground-state ansatz 



GS{ 



2pol.\ _ 



It) 



p\ 



+/r)+P2|+/r"-) -li) Pi -/r)+P2|-/r"-) 



(19) 



where |±/i) = e^-^^^^i |0). To optimise the state, we minimise the variational ground state energy E 



.2pol. 



Hi 



Qg2po\.\ ]2^j^gj.i(3g^}}y^ subject to the normalisation constraint 2p1 + 2p2 + 4pip2e 2(/f 



The Silbey-Harris state is again obtained by letting p2 

1 r, . 



0, such that 



= 



V2 



i+/r)-it)i-/r)] 



ranti . \ 2 

Jl ) = 1^ 



(20) 



and there is only a single displaced oscillator associated with each spin. Minimisation of E^^ = (^GSf^ \ Hi \GSf^)^ 

leads to a self-consistent equation for the optimised displacement, /f^ = g/[2{Ae~'^''^^^^^ +^i)]. The non-optimal 
bare polaron state has the same form as Eq. (20), but with the displacement fixed at fi = g/(2uoi). 

In Fig. [l] we plot the dimensionless ground state energy E/uji determined from our polaron- antipolaron variational 
ansatz as a function of the dimensionless spin-oscillator coupling strength g/uoi^ and compare with the results from 
an exact numerical diagonalisation of the model, and from Silbey-Harris and polaron theories (Eqs. ([6| and ^ in the 
single mode case, respectively). In this figure, A/cJi > 1 for all plots, and so we would expect standard polaron theory 
to break down in this regime, since the full oscillator displacement is no longer appropriate. From the dashed curves, 
this can indeed be seen to be the case, and polaron theory may even predict the incorrect trend as g/uji increases. 
Silbey-Harris theory fixes this problem to a certain extent (at least at small g/uji^ see dashed-dotted curves), though 
again runs into problems as the coupling strength increases, deviating from the numerically-exact results, and even 
more worryingly predicting discontinuous behaviour in the ground-state energy at certain values of g/uoi. Our polaron- 
antipolaron variational ansatz, however, predicts ground-state energies in almost perfect agreement with the numerical 
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FIG. 1. Ground state energy of the Rabi model. The single- mode spin boson model is considered as a function of 
spin-oscillator coupling strength, and the ground state energy calculated from our variational two-polaron ansatz (red solid 
curves), the Silbey-Harris one-polaron variational state (blue dashed curves), and from the bare polaron ground state (green 
dashed curves). Exact energies calculated by numerically diagonalising the Hamiltonian (for a basis of 200 states) are shown 
as grey dots. The values of A/cJi used are shown on each plot. For all plots, A/cJi > 1, which is the regime in which we expect 
our variational state to outperform the Silbey-Harris or bare polaron treatments. 



results for all possible coupling strengths. Furthermore, the discontinuous behaviour seen in Silbey-Harris theory is 
removed in this more flexible variational state. 

The failure of the single-polaron theories can be interpreted by analysing, in position space, the oscillator states 
associated with the spin projetions |^) and |t) in the model ground state {(l)^{x) and 0^(x), respectively), where 

we take the ground state to have the form \GSi) = |^) + |^), i.e. we absorb any normalisation factors 



and minus signs into the oscillator states, 



is thus negative using this convention. Shown in Fig. are thus 

0j(x) (red) and (/)^(x) (blue) as a function of position, calculated from our polaron- antipolaron ground state (solid 
curves), from the Silbey-Harris ground state (dashed-dotted curves), and from a numerical diagonalisation of the full 
Hamiltonian (points). We take A/cJi = 4 here as a representative example. For g/uji ~ 1, the displacement from 
X = is found to be fairly small, |/| < g/(2uji)^ which is why the full polaron approach fails, see Fig. [!}. The 
correct displacements can be captured by the Silbey-Harris theory and as well by the more flexible two-polaron ansatz 
presented in Eq. 



nature of the osci 



(19). However, as the coupling strength increases further, we can see that the double displacement 
lator states starts to become extremely important. For example, at g/uoi = 3 we observe that the 
Silbey-Harris state is completely unable to reproduce the correct oscillator wavefunctions, due to the restriction to 
a single displacement associated with each spin state. In fact, for these parameters, the displacements obtained by 
the Silbey-Harris approach are much too small, and reproduce none of displacements seen in our polaron-antipolaron 
ansatz, which itself matches the numerical solution very well. Finally, as the coupling strength is increased further, 
the Silbey-Harris displacements eventually "jump" to those of the full polaron transformation, which then captures 
the dominant displacements in the exact states quite well (see the g/uoi = 4 plot ), but of course completely misses 
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FIG. 2. Ground state wavefunctions of the Rabi model. Oscillator states 4>i{x) (red) and (blue) plotted as a 

function of position for A/cJi =4. Shown are the predictions of our variational ground state (solid curves), of the Silbey-Harris 
ground state (dashed-dotted curves), and exact calculations from a numerical diagonalisation of the full Hamiltonian (points). 



the smaller displacements in the opposite direction, which are still captured extremely well by our ansatz. Hence, the 
obtained ground state energy is still lower in our ansatz (and in the numerical diagonalisation) than from the Silbey- 
Harris state. The theories will eventually converge at larger coupling, however, when associating a single (polaron) 
displacement with each spin state finally becomes a good description. 



III. TWO-MODE RABI MODEL: ENTANGLEMENT ENTROPY 



In order to understand the non-monotoic behaviour of the joint spin-mode entropy 6'spin+/c results presented for 
the continuum environment in Fig. 5 of the main text, we consider in this appendix the entanglement between modes 
in the environment for the simpler case of a two-mode environment. The Hamiltonian for this case, the same used to 
generate the two- mode results in the main text, is given by 

H = ^(Jx - y XI ^^(^^ + + ^A^i' (21) 

i=l,2 i=l,2 

For this two- mode environment the ground state can be found by exact diagonalisation (ED) techniques, allowing a 
comparison to be made with the predictions of our anti-polaron ansatz. Figure [sjA shows results for 5'spin+i (we trace 
over mode 2) as a function of uj\ for ground state obtained by ED and three variational ansatze: Silbey-Harris (one- 
polaron), two-polaron, and three-polaron (to be discussed below) trial states. In all cases, a fixed relative detuning 
of the modes and ratio of coupling to frequency was kept, with CJ2 = l-OScJi, gi = 92 = 2.5cJi. For large frequencies 
uji ^ Ar, where we expect adiabatic polaron theory to describe the state, the product state (unentangled) form of the 
wave functions of oscillators 1 and 2 in the spin-projected states is expected to lead 6'spin+i to be controlled only by 
the polaronic correlations between mode 2 and the spin. This is determined by the renormalisation scale A^^, which is 
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a function of only the ratio 5^2 /^2- As this is held fixed in these simulations, we expect S'spin+i to approach a constant 
value at high frequency. This behaviour is indeed observed on all curves in Fig. [3]A., with the relative strong coupling 
{92/^2 — 2.5) leading to an almost fully- mixed spin state. According to Silbey-Harris theory, which has the same 
wave function structure as the adiabatic polaron theory but with different displacements, oscillators with frequencies 
well below A^^ should have suppressed displacements. Consequently, the renormalisation of the spin tunneling by slow 
modes should be continuously suppressed, leading to reduction of 5'spin+i as uj\ 0. This behaviour is precisely what 
is observed for the Silbey-Harris results in Fig. [3]A., with 6'spin+i decreasing monotonically with decreasing frequency 
(though with a sharp suppression below the spin-tunneling frequency scale). 




FIG. 3. Joint spin-mode entropies and oscillator wavefunctions for the two-mode spin-boson model. A. Joint 
spin-mode entropy 5'spin+i for a two-mode environment. For all these data (jJ2 = l.OScJi, gi = g2 = 2.5cJi and A = 0.01, 
sweeping the frequencies cji of the first mode. Plot shows results for (top to bottom) exact diagonalisation (grey stars), three- 
polaron ansatz (green diamonds), two-polaron ansatz (yellow squares) and one-polaron Silbey-Harris theory (red dots). The 
horizontal line indicates the maximum entropy of a fully-mixed spin state (ln(2)). The main panel is a close-up on the entropy 
peak occuring near the resonance frequency, as discussed in the text, while the inset shows the whole entropy and frequency 
range. B. Spin up-projected two- mode oscillator wave functions along the diagonal coordinate xi = X2 for the same parameters 
as the left panel, computed for cji = 0.0023 (namely at the peak position of the exact diagonalisation entropy curve in A). 
Results are shown for ground states obtained by exact diagonalisation (dashed black line), one-polaron Silbey-Harris ansatz 
(solid red line), two-polaron ansatz (solid yellow line) and three-polaron ansatz (solid green line). 



However, as in the multimode case, we see that the numerically exact ED ground state shows an entropy peak 
(in excess of the entropy of a fully-mixed spin state) at frequencies close to the scale A^^. As was shown in the 
two-mode wave function plots in the main text, the spin-projected oscillator wavefunctions in the antipolaron regime 
are entangled by their non- adiabatic responses, and these inter- mode correlations lead to the entropy peak seen in 
Fig. |3]A (due to both spin and mode 1 states becoming mixed when mode 2 is traced over). As in the single mode 
case, this occurs in the regime where antipolarons form, gi ^ coi ^ Aji. The absence of the antipolaron component in 
the Silbey-Harris theory means that this feature cannot be described by this state. 

However, a variationally-optimised two-polaron ground state ansatz as in Eq. ( 19) captures this peak structure in the 



frequency dependence of 5'spin+i- To understand qualitatively how antipolaron components in the two-polaron ansatz 
generate inter-mode entanglement and, thus, the increased spin-mode entropy, consider the (spin-up) spin-projected 
part of the ground state. This gives a contribution to the total density matrix of 



m,i,2 = I t)(t I (p?i/r'-)(/r'i +P2i/r"-)(/r"i +PiP2i/r'-)(/r"i +piP2i/r"-)(/r'i) 

where |/^) = e^'«=i'2 '^^^"^''^ |0). Tracing over mode 2, the reduced state p^^^i is 



m,i = I txt I (p?i/r'-)(/r'i +Pii/r"->{/rH +PiP2*22i/r'-){/f 



■PiP2*22i/r"->(/r' 



I) 



(22) 



(23) 



where ^22 = Tr2 



\f'. 



As /2 and /f"*'' have opposite signs when mode 2 is in the 
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antipolaron regime, the overlap integral ^22 suppresses the purity of the spin-projected state of mode 1. In the 
extreme case where ^22 ^ 0, the mixed reduced state is 

Pn,i ^ I t)(t I (p?i/r'-)(/r'i +P^i/r"-)(/r"i) • (24) 

This expression provides important intuition as to frequency dependence of inter-mode entanglement. The state 



given in Eq. (24) will have high entropy if the antipolaron weight P2/P1 in the wave function is significant and the states 
l/f"*^-) and l/p-) have weak overlap ( i.e. are close to orthogonal) . We therefore expect that the entropy peak will 
appear in the strong antipolaron regime (where P2/P1 and the displacements ~ _yanti. sizeable). However, 
this also requires that ^22 is also ~ 0, requiring that mode 2 is also in the strong anti-polaron regime. Our theory 
thus establishes the microscopic link between the appearance of intra-mode non-classicality (cat-states/antipolaron) 
and many-body entanglement between such modes. For the results of Fig. [sJA., the 5% detuning between modes 
leads to both modes showing antipolaron features at similar frequencies, leading to the large entropy peak. In the 
multimode case, we therefore expect to find antipolaron components in the environmental wavefunctions everywhere 
inside the region of positive 6'spin+/c — 5'spin, allowing the frequency range and position of the non-classical (cat-like) 
environmental features - which could be experimentally probed through the environmental response function - to be 
inferred. We have also checked that when the detuning between modes is made larger, such that the modes do not 
both develop antipolarons for the same parameters, the entropy peak becomes much smaller (not shown). 

Finally, we note that while the two-polaron ansatz captures the essential physics of the entropy peak, the preservation 
of (ax) at strong coupling, and also provides an excellent description of the shapes of the entangled wave functions, the 
agreement with the ED entropy is not perfect, with the two-polaron ansatz slightly underestimating the peak entropy, 
see Fig. [3]A.. To investigate this remaining small discrepancy, we have also implemented a variational three-polaron 
ansatz of the form: 

lG53p°i-) = 1 1 ) ® [pi|+/r') +P2|+/r"-) +P3I+/I)] (25) 
- 1 i ) ® [pi\-fr' )+P2\-fr'-)+P3\-f!)\ ■ 

Figure [3^ shows the wave functions of the spin up-projected states of the oscillators along the diagonal coordinate 
xi = X2 for the ground states obtained by ED and the two and three-polaron ansatze. Comparing the ED and 
two-polaron wavefunctions, we see that the two-polaron wavefunctions captures the displacements and weights of 
the polaron and antipolaron very well, but slightly underestimates the amplitude of the wavefunction around the 
origin. For simplicity, the three-polaron solution was determined by fixing /| = and treating all other parameters 
variationally. The result, shown in Fig. J3^ has almost perfect overlap with the ED results and gives an improved 
prediction for the entropy peak in Fig. [3]A.. This result suggests that it may be fruitful to consider a multipolaron 
expansion of the state in the many- mode cases, particularly if one is interested in reproducing sensitive measures of 
the many-body state structure (such as the joint entropy or other tomographic objects) rather than the simple spin 
observables (which are already well-approximated by the two-polaron results). 

IV. MULTI-MODE SPIN-BOSON MODEL: WIGNER DISTRIBUTIONS 

We discuss here various Wigner distributions associated to the reduced density matrix living in the subspace spanned 
by the qubit and one given oscillator mode with frequency uj^. The qubit degrees of freedom can be used for filtering 
out the polaron and antipolaron contributions within the wavefunction, thanks to appropriate insertions of Pauli 
matrices in the standard definition of the Wigner functiorPI. For instance, we can project onto the |t) component 
only, by considering: 

W[%{X) = e^(^-^)(G5|e^4-Aa,i + ^|^^\ (26) 

J TT^ 2 ' ' 

This Wigner distribution can be addressed from our two-polaron variational state In the case A <C cjc, we find 
P2 oc A/ujc ^ Pi^ and the wavefunction is normalized by taking simply pi 1/a/2. We thus need to compute: 

2 Z J TT I 

(27) 
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FIG. 4. Wigner distributions for the many-mode spin-boson model. A. Diagonal Wigner distribution defined in 



Eq. (26) computed by the NRG (dashed line) and compared to the one and two-polaron results found from Eq. (31 ). Parameters 
are a = 0.8 and A/cjc = 0.01. The one-polaron Silbey-Harris state is enough to account very well for the polaron content 
of the exact wavefunction. B. Off-diagonal Wigner distribution defined in Eq. (33) computed by the NRG (dashed line) and 
compared to the one and two-polaron results found from Eq. (34). The general magnitude and non-Gaussian form of the Wigner 



distribution is only accounted for by the two-polaron state, with a complete failure of the Silbey-Harris wavefunction (note the 
100 X magnification used to reveal its tiny contribution to the Wigner distribution). The remaining quantitative deviations in 
the two-polaron ansatz are due to the presence of additional antipolaronic contributions in the total wavefunction, as can be 
inferred from the computation done within a three-polaron trial state. 



Using usual coherent state algebra, we readily obtain the required overlaps: 
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and performing the Gaussian integral in Eq. ^71 yields 



w[%{X) = ^e 



-2(X-/P°^-)' 



P2g-2(X-/r*^-)' 
TT 



(31) 



Note that the second term in Eq. (31) provides only a small correction to the first contribution, a feature which 
follows from (i) p2 <^ 1 and (ii) the overlap appearing in this second contribution can be approximated as 
g-i Eg^fc(/r -/r'")^ - I /g - 6"^^'?^'^^''° '^^ Ai^/A < 1 (for a > 0.5) because the antipolaron is equal 

and opposite to the polaron at low energy. Similarly, the third term in Eq. ( [31] ), which would peak at the antipolaron 
displacement, is of order P2 <^ 1^ and so also provides a tiny contribution. Thus, the |t)-projected Wigner function 
is dominated by the purely polaronic contribution, as we indeed demonstrate by the impressive agreement with the 
numerically exact NRG computation of W^j+l^z (X) in the left panel of Fig. jij 

In order to highlight the emergence of antipolaronic contributions in the wavefunction, we now insert the off-diagonal 
(7+ Pauli matrix, which leads to the equivalent expression as defined in Eq. (7) of the main text: 



(32) 



We thus need to compute 
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A computation similar to performed above leads to the final result: 
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(34) 



The above expression shows important differences from the |t)-projected Wigner distribution of Eq. (31). Indeed, the 
first term associated with the purely polaronic response is now of order ^q^kUq ) ^ \— fq^^') — A^^/A <C 1, 

and so is subdominant to the second contribution (with mixed polaron-antipolaron origin) of order p2 oc A/lOc (Note 
that the overlap ^q^kifq +fq ) appearing in the second term is of order 1). The third term, of order P2A11/ A 
is even more smaller. Thus, the off-diagonal Wigner function W^^ {X) can be used to highlight the emergence of 
antipolarons in the many-body ground state wavefunction of the continuous spin-boson model, as was also discussed 
in the main text. 
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